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Abstract. A new Monte Carlo algorithm for 2-dimensional spin glasses is presented. The use of clusters 
makes possible global updates and leads to a gain in speed of several orders of magnitude. As an example, 
we study the 2-dimensional ±J Edwards- Anderson model. The new algorithm allows us to equilibrate 
systems of size 100 2 down to temperature T = 0.1. Our main result is that the correlation length diverges 
as an exponential (£ ~ e 2/3J ) and not as a power law as T —* T c = 0. 

PACS. 75.10.Nr Spin glass and other random models - 02.70.Lq Monte Carlo and statistical methods 



1 Introduction 

1.1 General considerations 

The understanding of (disordered) Ising ferromagnets has 
been greatly enhanced by fast Monte Carlo (MC) simula- 

' tions using cluster algorithms |Q,@|. Unfortunately those 
techniques cannot be directly applied to models such as 
spin glasses (SG) because of frustration. Attempts have 
been made to generalise cluster methods but the 

, resulting algorithms are complicated and the speed in- 
crease is not impressive. A cluster algorithm for fully frus- 

| trated systems already exists 0j but cannot tackle disor- 
der. Other techniques such as exchange MC (EMC) (also 
called parallel tempering) Q] allow big improvements over 

, standard one-spin flip MC and are widely used for SG. 
Nevertheless the sizes and temperatures accessible to sim- 
ulations are still not enough to clearly solve many im- 

■ portant issues (see || for a review). In the case of 2- 
dimensional spin glasses, the best method to date is the 
replica MC (RMC) by Swendsen and Wang || (which es- 
sentially reduces to EMC in higher dimensions). 

We present here a new cluster MC algorithm for 2- 
dimensional SG which is much faster than previous MC 
techniques (namely RMC). It thus gives access to sizes 
and temperatures which were unreachable before. Note 
that transfer matrix methods which are widely used for 
2-dimensional systems are limited to "small" sizes, usu- 
ally no more than 16 x oo which appears to be not enough 
to tackle certain open problems. With our new tool we 
have studied the SG transition in the 2-dimensional ±J 
Edwards- Anderson (EA) (To) model for which several ques- 
tions arc still unsettled. In particular the value of the criti- 
cal temperature (zero or not) |Tl], |I^,[l3,[l4, D| and the na- 
ture of the divergences (power laws or exponentials) jl6| 
are still debated. We present evidence at the end of this 
article that T c = and that the correlation length follows 



an exponential law £ ~ e 2 @ J , which is different from the 
standard lore. 



1.2 Models 

The model we consider here is a general Ising spin model 
on a lattice, whose Hamiltonian is: 



H{S) 



/ JjjSj Sj ^ ' hjSj, 



(1) 



where Si = ±1. The interactions Jy and the magnetic 
fields hi are any fixed real numbers. Let be the number 
of spins. Although the algorithm described below is cor- 
rect in all cases, it is faster than more traditional methods 
only in the 2-dimensional case with nearest neighbour in- 
teractions. This allow to study the Ising spin glass, the 
(disordered) Ising ferromagnet and the random field Ising 
model in two dimensions. We will only consider those cases 
in the following. 



2 The algorithm 

2.1 The cluster Monte Carlo step 

Let us now describe a Monte Carlo move which allows a 
global update of the spin configurations. Consider a sys- 
tem consisting of two independent spin configurations at 
the same temperature 1/(3; thus a "configuration" of this 
system is a set of two spin configurations: C = {{S}}, {Sf}). 
The same Hamiltonian (the same J^ 's and h^s) apply to 
both configurations. We want to sample the configurations 
with the weight: 



P{C) oc exp [-0 (HiS 1 ) + H(S 2 ))] 



(2) 
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Fig. 1. An example of the value of the q^s (1 is white, —1 is 
black). The system is here a 200 2 ±1 2-dimensional Edwards- 
Anderson spin glass at equilibrium at /3 — 1.75. 



To do this, it is sufficient to have an ergodic Markov chain 
and to enforce the "detailed balance" condition: 



P{C)n t 



P(C)n t 



C'->C, 



(3) 



where Uq^c is the transition probability from C to C. 

We define the local overlap at site i between the two 
replicas by = SfSf. This defines two domains on the 
lattice, the sites with % = 1 and the sites with qi = — 1. 
We call clusters the connected parts of these domains (two 
sites i and j are connected if Jij ^ 0) . Our cluster Monte 
Carlo step proceeds as follow: choose one site at random 
for which q{ — — 1, and flip the cluster to which it be- 
longs in both configurations, ft is quite easy to check that 
iJ(S' 1 ) + H(S 2 ) is unchanged by this transformation, the 
interface and magnetic energies of the cluster in both con- 
figurations being exchanged by the flip. Moreover the (ft's 
are also unchanged, together with the definition of the 
clusters. It is thus clear that equation || is verified for this 
step since P(C) = P(C) and 7T c ^c = Ptc^c- Note that 
in the case where the /i,'s are zero, there is no magnetic 
energy and one can also choose a spin with q L = +1 to flip 
a cluster. In figure f] an example of the values of the q^s 
is presented. The clusters are the connected parts of the 
white and black domains. 

Our construction is very similar to RMC M with the 
essential difference that both replicas have the same tem- 
perature. This point is very important for two reasons: (i) 
our cluster moves are accepted with probability one, (ii) 
we can use more than two configurations as explained in 
the next section. The definition of the clusters is also rem- 
iniscent of |l7| where a cluster a la Wolff is grown inside a 
region where two replicas are different. This last method 
works in random field but requires non frustrated inter- 



actions (as far as we know it has never been extended to 
frustrated interactions). 



2.2 Description of the algorithm 

In spin systems one is essentially interested in three ob- 
servables: the energy of the system (H), the magneti- 
sation (m = Y^Si/N) and for spin glasses the overlap 
(q = J2SfSi/N)- But from the previous discussion it is 
clear that the proposed cluster flip keeps all those quan- 
tities constant for the whole system though the energies 
and magnetisations of each replica do change. So as de- 
scribed in the previous section the algorithm is non er- 
godic. To build a valid Monte Carlo algorithm, we need 
first to enforce ergodicity. The simplest way to achieve this 
is to add another kind of move: We choose the standard 
one-spin flip move with Metropolis acceptance (as done in 
RMC). The main point of this new algorithm is that we 
do not restrict ourselves to only two configurations. We 
use n>2 replicas at the same temperature. This trick al- 
lows a much faster relaxation because these n replicas are 
mixed together very quickly. In particular with 3 replicas, 
if one flips a cluster between replicas 1 and 2, even if q±2 
is unchanged, qis and 923 do change. 

At this point one problem remains, namely the total 
energy of the n replicas is conserved by the cluster moves. 
Only the one-spin flip moves can relax the total energy 
and the resulting algorithm is much slower than RMC 
(and also slower than EMC). That is why we embed our 
cluster method in an EMC. 

Here are the details of our algorithm: the system is 
composed of m sets at different temperatures. Each set 
consists in n independent replicas at the same temper- 
ature. The same Hamiltonian (with the same Jy 's and 
hi's) applies to all configurations. The algorithm repeat- 
edly does the following: 

1. One-spin flips: Do a standard one-spin flip move for 
each spin in each replica. 

2. Cluster moves: For each temperature, randomly par- 
tition the n replicas in pairs and do one cluster move 
for each pair. 

3. EMC: For each pair of neighbouring temperatures, do 
n standard EMC updates between the two sets (pairing 
each replica from one set with one from the other). 

In a standard EMC update, two spin configurations at 
different temperatures are exchanged with probability 



Pu 



min ( 1 , e 



=,(/32-/3i)(if2-fli) 



(4) 



The choice of the m temperatures is made as in EMC: 
One needs enough temperatures for the energy distribu- 
tions of two neighbouring temperatures to overlap; then 
the exchange moves can be accepted with a reasonable 
probability. A good rule of thumb is to try to fix the ex- 
change acceptance rate to 1/3 to optimise the mixing ef- 
fect. The choice of n seems to be simple: the larger the 
better. In fact, the computation time increases linearly 
with n but one obtains n independent configurations at 
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somehow the EMC is able to explore the energy landscape 
vertically (up and down in energy) whereas the clusters 
allow an horizontal search (search at constant energy). 
The cluster moves give a very quick mixing of the differ- 
ent configurations at one temperature and EMC allows 
for a renewal of the population. Even if the clusters are 
the same as in RMC, RMC does not perform so well be- 
cause of a less good mixing: In RMC the cluster moves 
always happen between the same configurations and they 
are hampered by the difference in temperature. 
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2.3 Only in 2 dimensions 



Fig. 2. Relaxation of the energy in a 100 2 spin glass for dif- 
ferent algorithms at /3 = 10. See the text for details. 
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Consider what happens in 3 dimensions: During the clus- 
ter construction, usually there are about as many sites 
with qi — 1 as with qi = — 1. In d = 3 the site percolation 
threshold is roughly 0.3 and so both q = 1 and q = — 1 sites 
percolate forming two system-size clusters (and some very 
small clusters). It is easy to see that nipping one of those 
big clusters in both replicas is the same as exchanging the 
replicas (except for the very small clusters). Hence, the 
move does essentially nothing (for the same reason, RMC 
becomes equivalent to EMC in 3 dimensions). This prob- 
lem is encountered as soon as the site percolation thresh- 
old is less than 0.5 which essentially forbids all lattices 
of dimension d > 2 or with more than nearest neighbour 
connections in d = 2. 



3 The EA model in d 



Fig. 3. Relaxation of q 2 in a 100 2 spin glass for different algo- 
rithms at P = 10. See the text for details. 

each step for each temperature. On the other hand the 
mixing effect of the clusters rapidly grows with n, so a 
large value is preferable, In the following we set n = 32, 
the maximum value allowed by our implementation (spin 
coding on 4- byte integers); it would probably be better to 
use even larger values of n. 

We investigate the relaxation time on a 100 2 , ±1 spin 
glass (hi — 0) with periodic boundary conditions at (3 — 
10 using n — 32 and m = 26 temperatures in the range 
(3 = 1 ... 10. We compare our algorithm to EMC and RMC 
(using the same temperatures) starting from random con- 
figurations and averaging over 32 independent runs. In 
figures || and ||, one can see the dramatic improvement 
over the EMC and the big improvement over RMC (at 
least a factor 100 in this case). Note that to determine if 
equilibrium has been reached, we also reran the algorithm 
starting from ground state configurations (obtained from 
an improved version of the algorithm described in |lS}| ). 
The equilibrium value is reached where the curves meet. 
It appears that the gain in speed increases with N and (3 
and may become really huge. 

Our algorithm is very efficient because the two tech- 
niques (clusters and EMC) are perfectly complementary: 



Using our algorithm, we have studied the 2-dimensional 
J = ±1 EA model with periodic boundary conditions on 
square L x L lattices. We have simulated 4 sizes (L = 
12, 25, 50, 100) for temperatures ranging from (3 = 0.3 
to j3 = 10. We have performed a disorder average over 
many samples (respectively 400, 400, 200 and 100). The 
large sizes and low temperatures involved made it difficult 
to equilibrate the system so we proceeded as follows. For 
each sample we first started from random configurations 
and waited for the lowest temperature to be completely 
occupied by ground states. Then we used these ground 
states as initial configurations for all temperatures and 
let the algorithm run. The point here is that it is easier 
for the algorithm to heat the system than to cool it, since 
the main obstacle is the low entropy of low energy states. 
We respectively used a relaxation time of 10 3 , 5.10 3 , 2.10 4 
and 10 5 MC steps for the different sizes. 

In each case we gathered statistics for the overlap q 
and measured the SG susceptibility 

X - NW), (5) 

where the over-line denotes the disorder average and the 
brackets the thermal average. We also measured the Binder 
cumulant 

'-K*-w> 
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Fig. 5. The SG susceptibility \ as a function of /3. Error bars 
Fig. 4. The Binder cumulant g as a function of the inverse are smaller than the symbols, 
temperature /3. The error bars are shown only for L = 100, 

they are smaller than the symbols for the others sizes. 1 i 1 1 1 1 ■ 



The standard expected scaling forms JllJ] are £ ~ (7 1 - 
T c )~ v for the correlation length, 



o-a L = 


12 


D — □ L = 


25 


L = 


50 


nL = 


100 



X^L 2 ^x(L l/u {T-T c )) 
for the susceptibility and 

g~g(l}/ v {T-T c )) 



(7) 

(8) 



for the Binder cumulant. Then the value of g should be 
constant at the transition and the curves in figure | should 
cross at T = T c . The curves definitely cross, but the cross- 
ing point goes to higher and higher values of (3 as the size 
increases which is a strong indication that no finite tem- 
perature transition occurs (i.e. T c = 0). Even if all the 
curves have a similar shape, it is interesting to note that 
the plateaus at low temperature cannot simply scale as 
predicted by equation || (nor by equation |l0|) , since those 
equations predict that g(L,T = 0) should be a constant 
(for T c = 0) which is clearly not the case. This behaviour 
is probably a finite size effect due to the fact that at these 
temperatures the system is in its ground state. The results 
for x ar e presented in figure ^[ 

Figures Band [?] show the result of the scaling of equa- 
tions and]| using \ jv = 0.35 and rj — 0.2 (and T c = 0). 
In both figures one can see a clear trend with size and 
the curves do not to overlap. This effect is not due to a 
bad choice of the parameters since the curves on figure [3] 
cross and no horizontal scaling around could make the 
curves overlap. The same argument applies to figure [?| for 
both axes. For the same reason corrections to scaling as 
proposed in JIJ] cannot improve the scaling. This leads us 
to conclude that this scaling does not apply. 

A few years ago, Saul and Kardar (TJ| proposed a few 
years ago an exponential scaling for the correlation length, 
namely £ ~ e 2/3 , which leads to: 




Fig. 6. The power law scaling: g as a function of TL 
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(9) 



Fig. 7. The power law scaling: x/L 1 ' as a function of TL 
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Fig. 8. The exponential scaling: g as a function of /3 — i \nL 
(no free parameters). The insert shows an enlargement of the 
high temperature region. 
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4 Conclusion 

We have described a cluster Monte Carlo algorithm for 
Ising spin models which incorporates an exchange Monte 
Carlo along with constant energy "cluster" moves. In 2 
dimensions with nearest neighbour interactions, this algo- 
rithm allows a gain in speed of several order of magnitude 
as compared to replica Monte Carlo and exchange Monte 
Carlo. And it works for any kind of interactions and mag- 
netic field (even frustrated and disordered). 

The ±J 2-dimensional Edwards- Anderson Ising spin 
glass was then studied at very low temperatures and for 
large sizes. The critical temperature appears to be exactly 
T c — and the standard power law scaling must be re- 
placed by the exponential sca ling proposed by Saul and 
Kardar (equations || and |ic| ). This probably means 
that the lower critical dimension is di = 2 for spin glasses. 
In an early work JljJ, McMillan stated that di > 2 and 

that at d = di one should have £ ~ e K @ ; we also tried this 
scaling law but it definitely does not apply in our case. 

As argued in (2(J excitations with non-trivial topology 
are necessary to have a spin glass phase, and it is pre- 
cisely in this case that the algorithm proposed here does 
not work. This point of view is also compatible with [ fl2"| 
where a 2-dimensional spin glass with next-nearest neigh- 
bour ferromagnetic interactions seems to have a finite crit- 
ical temperature (the percolation threshold is then less 
than 1/2 and excitations with non-trivial topology are 
possible). 

Only the ± J spin glass has been studied here; it would 
be interesting to see if the spin glass with Gaussian cou- 
plings behaves in the same way; this will be the subject 
of a publication to come. 
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Fig. 9. The exponential scaling: x/L 1 ' 8 as a function of 
P — \ InL. The insert shows an enlargement of the high tem- 
perature region. 

and 

ff ~5(/3-ilnL). (10) 

The corresponding scalings are shown in figures | and | 
with r\ = 0.2. The overlap is much better and there are 
fewer free parameters, which seems to indicate that this 
scaling is the correct one. In the inserts of these figures 
on can see that the discrepancies disappear as the size in- 
creases. Finally the scalings with T c ^ have also been 
tried (data not shown), but they give less good results 
than the exponential scaling though having two parame- 
ters more. 
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